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Abstract 


This  report  introduces  a  mixture  distribution  approach  to  modeling  the 
probability  density  function  for  lead  time  demand  (LTD)  in  problems  where  a 
continuous  review  inventory  system  is  implemented.  The  method  differs  from  the 
typical  “moment-matching”  approach  by  focusing  on  building  up  an  accurate,  closed- 
form  approximation  to  the  LTD  distribution  from  its  components  by  using  mixtures  of 
truncated  exponential  (MTE)  functions.  First,  construction  of  the  lead  time 
distribution  is  illustrated  and  the  approach  is  compared  to  two  other  possible  lead 
time  distributions.  This  distribution  is  then  utilized  to  determine  optimal  order 
policies  in  cases  where  a  buyer  makes  its  decisions  alone,  and  later  in  a  situation 
where  members  of  a  two-level  supply  chain  coordinate  their  actions.  Next,  a  mixture 
of  polynomials  (MOP)  approach  is  introduced  and  utilized  to  model  the  LTD 
distribution  for  the  case  of  discrete  lead  time  and  a  daily  demand  distribution  that 
assumes  a  standard  form.  Finally,  the  MOP  model  is  extended  to  model  an  LTD 
distribution  from  empirical  lead  time  and  daily  demand  data. 
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Modeling  Uncertainty  in  Miiitary  Suppiy 
Chain  Management  Decisions 


Numerous  probability  models  have  been  suggested  for  representing  uncertain 
demand  during  lead  time  in  continuous-review  inventory  management  systems  when 
both  lead  time  and  demand  per  unit  time  are  variable.  A  common  approach  to 
finding  a  distribution  for  lead  time  demand  involves  modeling  lead  time  (LT)  and 
demand  per  unit  time  (DPUT)  with  standard  probability  density  functions  (PDFs). 
Based  on  the  distributions  assigned,  a  compound  probability  distribution  is 
determined  for  demand  during  lead  time,  or  lead  time  demand  (LTD).  The  latter 
distribution  is  used  to  determine  reorder  point  and  safety  stock  policies,  and  may  be 
used  to  estimate  inventory  costs.  In  some  cases,  analytical  formulas  for  optimal 
reorder  point,  safety  stock,  or  stockout  costs  are  available  in  terms  of  the  compound 
distribution's  parameters,  while  in  other  situations  the  values  associated  with  certain 
percentiles  of  the  compound  LTD  distribution  are  estimated  to  provide  these  values. 
While  the  problem  of  finding  an  appropriate  LTD  distribution  has  been  well-studied, 
papers  written  in  recent  years  have  continued  to  pursue  methods  that  overcome 
unrealistic  distributional  assumptions  (Ruiz-Torres  and  Mahmoodi,  2010;  Vernimmen 
et  al,  2008). 

This  paper  illustrates  an  approach  for  constructing  a  mixture  distribution  for 
LTD  that  allows  the  LT  and  DPUT  distributions  to  be  state-dependent.  This  method 
also  allows  input  distributions  that  take  any  standard  or  empirical  form.  Use  of  the 
mixture  distribution  technique  is  first  demonstrated  in  the  context  described  by  Cobb 
(2013),  which  is  a  single  item  continuous-review  inventory  model  for  one  buyer.  For 
single-firm  operating  in  a  continuous-review  inventory  system,  the  mixture 
distribution  method  for  modeling  the  LTD  distribuition  differs  from  the  typical 
“moment-matching”  approach.  The  method  focuses  on  building  up  an  accurate, 
closed-form  approximation  to  the  LTD  distribution  from  its  components  by  using 
mixtures  of  truncated  exponential  (MTE)  functions. 

After  the  mixture  distribution  approach  is  described,  a  two-level  supply  chain 
model  where  the  buyer  operates  under  uncertain  demand  and  utilizes  a  continuous 
review  inventory  system  will  be  considered.  In  this  two-echelon  supply  chain  model, 
credit  terms  (Chaharsooghi  and  Fleydari,  2010),  quantity  discounts  (Li  and  Liu,  2006; 
Chaharsooghi  et  al.,  2011),  and  rebates  (Cobb  and  Johnson,  2014)  have  been 
suggested  as  coordinating  incentives  that  allow  the  supply  chain  members  to  divide 
the  cost  savings  resulting  from  coordinating  their  order  quantity  and  reorder  point 
decisions.  In  each  of  these  cases,  LTD  is  assumed  to  be  normally  distributed.  This 
assumption  is  not  always  realistic,  particularly  when  demand  per  unit  time  and  lead 
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time  are  each  random  variables  such  that  LTD  has  a  compound  probability 
distribution  (Eppen  and  Martin,  1988;  Lau  and  Lau,  2003;  Lin,  2008).  This  paper  will 
incorporate  the  previously  described  model  (Cobb,  2013)  into  the  two-echelon 
supply  chain  problem  to  show  that  this  model  can  obviate  the  need  to  assume  that 
demand  for  the  entire  lead  time  period  is  normally  distributed. 

Mixture  of  polynomials  (MOP)  models  (Shenoy  and  West,  2011)  are  an 
alternative  to  the  MTE  model  for  approximating  PDFs.  These  models  are 
implemented  in  the  two-level  supply  chain  model  in  two  situations.  First,  MOPs  are 
used  to  fit  LTD  distributions  given  each  possible  lead  time  value  when  the  PDFs 
have  a  standard  functional  form.  Next,  these  distributions  are  approximated  from 
historical  data.  In  each  case,  the  mixture  distribution  approach  can  be  applied  to 
calculate  a  closed-form  approximation  to  the  LTD  distribution. 

The  next  section  describes  lead  time  demand  distributions  and  uses  an 
example  dataset  to  show  how  standard  PDFs  can  be  used  as  approximations  to  the 
LTD  distribution.  The  mixture  distribution  method  is  also  used  for  the  example 
problem.  Next,  the  different  approximations  to  the  LTD  distribution  are  used  to  find 
optimal  inventory  order  quantity  and  reorder  point  policies.  This  is  followed  by  an 
illustration  of  how  the  mixture  distribution  approach  can  allow  more  complicated  LTD 
distributions  to  be  incorporated  into  such  problems.  The  two-level  supply  chain 
model  is  then  introduced,  and  the  mixture  distribution  approach  is  used  to  model 
LTD  in  the  context  of  decentralized,  centralized,  and  coordinated  supply  chains.  In 
the  next  two  sections,  the  MOP  approximations  are  described  for  the  standard  PDF 
case  and  the  situation  where  the  MOP  distributions  are  estimated  from  historical 
data.  The  final  section  concludes  the  paper. 

Lead  Time  Demand  Distributions 

LTD  in  a  continuous-review  inventory  system  is  often  assumed  to  follow  a 
compound  probability  distribution.  Suppose  L  is  a  random  variable  for  lead  time  (LT) 
and  D  represents  random  demand  per  unit  of  time  (DPUT).  LTD  is  a  random 
variable  ^determined  as 


X  —  D2  Di  Di  .  (1 ) 

Therefore,  Xis  a  sum  of  random,  independent  and  identically  distributed 
(i.i.d.)  instances  of  demand.  The  mean  and  variance  of  Xcan  be  calculated  as 

EiX)  =  E(_L)  ■  EiD)  and  Var(X)  =  E(_L)  ■  VariD)  +  [EiD)]^  ■  Var(L)  (2) 

Suppose  the  data  in  Table  1  is  available  to  estimate  a  LTD  distribution.  This 
table  contains  50  observations  of  daily  demand  for  an  inventory  item  and  10 
observations  for  lead  time  on  orders  of  the  same  item.  The  expected  value  of  daily 
demand  is  E(D)=2.88  and  the  variance  of  this  random  variable  is  Var{D)=2.84.  Lead 
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time  has  an  expected  value  and  variance  of  E{L)=5.3  and  Var{L)=6.9,  respectively. 
According  to  the  formulas  in  (2),  the  expected  value  and  variance  of  LTD  are 
E(X)=15.26  and  Var{X)=72.3,  respectively. 

The  remainder  of  this  section  will  illustrate  three  possible  methods  for 
approximating  the  LTD  distribution  underlying  the  data  in  Table  1 . 

Normal  Approximation 

The  service  level  is  defined  as  the  percentage  of  replenishment  order  cycles 
where  demand  during  lead  time  is  satisfied.  To  determine  the  reorder  point  {R) 
required  to  achieve  a  desired  service  level,  a  typical  textbook  approach  is  to  assume 
the  LTD  distribution  is  normal  and  use  normal  distribution  tables  or  Excel  formulas. 
For  example,  to  find  the  R  needed  to  achieve  a  95%  service  level  for  the  LTD 
distribution  with  expected  value  and  variance  described  in  Table  1 ,  the  Excel  formula 
NORM.INV(0.95,15.26,72.3''0.5)  can  be  used  to  find  F?=29.25. 


Table  1.  Observations  for  Daily  Demand  and  Lead  Time 


Daily  demand  (DPUT) 

1 

2 

2 

1 

4 

1 

1 

1 

1 

1 

3 

5 

3 

2 

5 

4 

2 

2 

3 

2 

2 

3 

3 

3 

1 

3 

6 

3 

6 

2 

5 

1 

5 

3 

2 

6 

1 

2 

4 

1 

3 

2 

2 

2 

6 

5 

5 

1 

3 

7 

Lead  time  (LT) 

3 

5 

3 

4 

4 

5 

5 

10 

5 

10 

The  normal  approximation  to  the  LTD  distribution  and  the  reorder  point 
R=29.25  are  illustrated  graphically  in  Figure  1 .  By  implementing  this  policy,  we 
would  expect  to  stockout  on  5%  of  replenishment  order  cycles. 
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Figure  1.  LTD  Distribution  and  Reorder  Point. 


Negative  Binomial  Approximation 

While  the  normal  approximation  to  the  LTD  distribution  is  popular,  there  are 
numerous  other  approximations  that  have  been  suggested  in  the  literature.  For 
example,  Taylor  (1961)  suggests  using  the  negative  binomial  (NB)  distribution  for 
the  case  where  the  Poisson  distribution  is  a  good  fit  for  DPUT  and  LT  has  a  gamma 
distribution.  Denote  the  approximate  LTD  distribution  by  /.  Here  we  assume  the 
NB(r,p)  distribution  for  LTD  is 

/(x;r,p)  x  =  0,1,2,...  (3) 

where  r(-)  is  the  gamma  function.  Given  this  formulation,  E(X)  =  rp/(l  -  p)  and 
Var(X)  =  E(X)/(1  -  p).  There  are  two  ways  of  finding  a  reorder  point  that  will 
provide  an  appropriate  service  level  with  this  NB  formulation.  Taylor  (1961)  provides 
a  formula  to  calculate  stockout  probabilities  as  a  function  of  the  underlying  Poisson 
and  Gamma  distributions.  These  can  be  calculated  for  possible  reorder  point  values 
until  a  suitable  value  that  meets  the  service  level  objective  is  found.  Excel  can  also 
be  used  to  enumerate  the  probabilities  of  achieving  a  certain  service  level  with 
various  possible  values  of  R.  Unfortunately,  the  built-in  NEGBINOM.DIST  function 
only  accepts  integer  values  of  the  r  parameter,  so  these  probabilities  must  be 
calculated  using  the  formula  in  (3)  and  the  GAMMALN  function. 

For  the  data  in  Table  1 ,  we  can  use  the  empirical  expected  value  and 
variance  to  solve  two  equations  and  two  unknowns  and  obtain  r=4.08  and  p=0.79. 
This  NB  distribution  is  shown  in  Figure  2.  The  value  of  F?that  provides 
approximately  a  95%  service  level  is  R=3^ . 
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Figure  2.  Negative  Binomiai  Distribution  for  LTD. 


This  solution  is  essentially  the  same  as  the  one  found  using  Taylor’s  (1961) 
analytical  formulas.  In  this  case,  the  Poisson  daily  demand  assumption  may  be 
reasonable,  because  E(D)  and  Var{D)  are  very  similar,  a  feature  of  the  Poisson 
distribution. 


Mixtures  of  Truncated  Exponentials  (MTE)  Approximation 

The  functional  form  of  some  PDFs,  such  as  the  negative  binomial  PDF  in  (3) 
do  not  permit  integration  in  closed-form.  The  means  the  result  of  an  expected  value 
calculation  with  such  a  PDF  does  not  have  a  functional  form  that  can  be  used  for 
further  computation.  These  calculations  could  include,  for  example,  building  a  cost 
function  to  perform  nonlinear  optimization  to  find  optimal  inventory  policies.  One 
approach  suggested  to  overcome  this  limitation  is  the  MTE  model  (Moral  et  al., 
2001). 

An  example  of  a  4-piece,  2-term  (ignoring  the  constant)  MTE  function  that  can 
be  used  to  model  LTD  given  a  lead  time  of  L=3  for  the  problem  in  the  previous 
section  is: 


/x|{L=3}W  — 


r-0.7148  +  0.6681  exp{0.0325x}  +  0.000048  exp{0.989x} 

I  -96.721  -  318.54  exp{-1.945x}  +  96.76  exp{0.000128x} 
I  0.1383  -  (1.63E  -  06)  exp{x}  +  (2.89E  -  09)exp{1.5x} 
1-0.0252  +  0.9786  exp{-0.205x} 


if  2.5  <  X  <  5 
if  5  <  X  <  8 
if8  <  X  <  11.5 
if  11.5  <  X  <  17.5 


(3) 


This  function  was  found  by  simulating  500  series  of  three  observations  for 
daily  demand  from  values  in  Table  1  using  a  bootstrapping  approach.  The 
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constants — coefficients  on  the  exponential  terms  and  coefficients  on  the  variable 
X— were  determined  by  fitting  a  function  to  the  simulated  histogram.  There  is  an 
established  literature  on  fitting  MTE  functions  to  historical  data;  in  this  case,  the 
method  suggested  by  Moral  et  al.  (2002)  was  utilized.  A  graphical  view  of  the  MTE 
function  overlaid  on  the  simulated  histogram  is  shown  in  Figure  3. 


Figure  3.  MTE  Distribution  for  LTD  Given  a  Lead  Time  of  3  days. 

Similar  functions  fx\{L=i}  can  be  constructed  for  the  other  possible  lead  time 
values,  L=4,  5,  and  1 0.  From  the  data  on  lead  time  observations  in  Table  1 ,  we  can 
estimate  P{L=3)=  P{L=4)  =  P{L=^0)  =  0.2  and  P(/.=5)=0.4.  A  mixture  distribution 
approach  (Cobb,  2013)  can  be  employed  to  find  the  LTD  distribution.  Flere,  the  LTD 
distribution  is  determined  as 

fxix)  =  P(L  =  3)  ■  fx\{L=3](.^)  +  —  4)  ■  fx\[L=4}(.^)  +  =  5)  ■  fx\[L=5](.^)  (4) 

+P(.L  =  10)  ■  fx\{L=10](.^)- 

The  MTE  function  is  shown  in  Figure  4  overlaid  on  the  previously  described 
NB  distribution.  This  MTE  function  has  17  pieces  and  up  to  six  terms  in  each  piece. 
For  illustrative  purposes,  a  continuous  NB  parameterization  is  displayed.  Since  the 
class  of  MTE  functions  is  closed  under  addition,  multiplication,  and  integration  (Moral 
et  al.,  2001),  the  mixture  distribution  resulting  from  the  calculation  above  is  also  an 
MTE  function.  Thus,  it  retains  the  same  desirable  mathematical  properties. 

We  can  perform  closed-form  integrations  of  the  MTE  LTD  distribution  to  find  a 
reorder  point  that  achieves  a  desired  service  level.  In  this  case, 

33.3 

j  /xW  ~  0-95  , 

0 

SO  we  can  set  R=33.3  to  obtain  a  95%  service  level. 
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Figure  4.  MTE  LTD  Distribution  Overiaid  on  a  NB  Approximation. 

The  next  section  discusses  the  use  of  the  MTE  function  for  finding  inventory 
policies  in  a  continuous-review  inventory  system. 

Calculating  Inventory  Policies 

Suppose  that  we  want  to  determine  an  optimal  order  quantity  and  reorder 
point  in  a  continuous-review  inventory  system  (a  “{Q,R)”  policy).  We  will  consider 
four  models  that  could  be  used  to  find  the  best  policy  given  the  data  available  (see 
Table  1):  1)  a  normal  approximation  to  the  LTD  distribution;  2)  the  NB  approximation 
to  the  LTD  distribution;  3)  the  MTE  mixture  distribution;  and  4)  a  simulation- 
optimization  model  that  simulates  lead  time  and  demand  values  from  the  empirical 
distributions  developed  from  Table  1 .  We  term  the  latter  model  the  “actual”  solution. 

A  simple  cost  function  with  no  backordering  allowed  (Johnson  and 
Montgomery,  1974)  for  this  problem  is 

TC(Q,R)=K-^  +  ^^  +  h-(O.SQ+R-E(X))  .  (5) 

In  this  equation,  K\s  the  fixed  cost  per  order,  Y\s  the  expected  annual  demand,  h  is 
the  holding  cost  per  unit  per  year,  and  n  is  the  stockout  cost  per  unit.  The  average 
inventory  includes  safety  stock  of  R-E{X).  The  shape  of  the  distribution  for  LTD 
determines  the  expected  shortage  per  cycle,  Sr.  For  a  given  reorder  point, 

Sr  =  J“(x  -  R)  ■  fxix)  dx.  (6) 

Suppose  Y=E{D)  ■  250  working  days  =  720,  K=30,  h=4,  and  tt=5.  The  key  to  finding 
an  optimal  {Q,R)  combination  is  to  evaluate  Sr  as  part  of  constructing  the  total  cost 
function  in  (5).  With  the  MTE  function,  the  calculation  in  (6)  can  be  performed  in 
closed-form,  and  the  result  substituted  into  (5)  to  obtain  a  closed-form  total  cost 


.miiWTiArniQifiim, 
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function.  The  expected  shortage  per  cycle  as  a  function  of  R  is  an  8-piece 
expression,  with  selected  terms  shown  below: 

SR(r) 

^-3876.5  +  4.66  exp{-0.205r}  +  6.31  exp{-0.172r} 

+3888.1  exp{0.005r}  -  21.82r  -  0.04r^  if  16.15  <  r  <  16.5 

-3890.6  +  4.66  exp{-0.205r}  +  6.31  exp{-0.172r} 

+3888.1  exp{0.005r}  +  20.6  exp{-0.140r}  -  20.64r  -  O.OVr^  if  16.5  <  r  <  17.5 
'  -3889.2  +  6.31  exp{-0.172r}  +  3888.1  exp{0.005r} 

+20.6  exp{-0.140r}  -  20.76r  -  0.07r^  if  17.5  <  r  <  23.5 

.-8.74  +  29.87exp{-0.78r}  +  0.28r  -  0.002r^  if  31  <  r  <  46.5  . 

Optimization  over  the  resulting  cost  function  is  fast.  The  example  here  was  solved  in 
Mathematica  9.0  by  using  the  ArgMin  function.  The  results  obtained  using  the  four 
methods  under  consideration  are  shown  in  Table  2.  An  iterative  approach  (Hadley 
and  Whitin,  1961)  in  combination  with  numerical  integration  was  implemented  to  find 
the  solutions  using  the  normal  or  NB  approximations.  The  table  shows  the  values 
Q*  and  R*  which — when  implemented  simultaneously — minimize  annual  total  cost. 
The  computing  (CPU)  times  required  to  obtain  the  solutions  are  also  shown.  The 
simulation-optimization  solution  was  simply  stopped  after  running  for  several  hours, 
and  the  values  obtained  were  assumed  to  be  the  best  possible  solution. 

Table  2.  Results  for  Inventory  Policies  Determined  Using  Four 


Approaches. 


Method 

Q* 

FT 

TC 

CPU  (sec.) 

Normal  Approximation 

108 

25 

482.99 

3.57 

NB  Approximation 

110 

25 

482.89 

3.76 

MTE  Mixture  Distribution 

110 

27 

481.10 

1.26 

Simulation-Optimization 

108 

27 

480.82 

oo 

Table  2  shows  that  the  MTE  mixture  distribution  works  equally  as  well  as  the  other 
approaches  when  implemented  to  obtain  an  optimal  {Q,R)  policy.  The  next  section 
illustrates  that  the  mixture  distribution  approach  can  be  used  to  model  more 
complicated  LTD  distributions. 

State-Dependent  Variables 

The  advantage  of  the  mixture  distribution  approach  (Cobb,  2013)  in  inventory 
management  problems  is  that  more  complex  LTD  distributions  can  be  constructed 
by  building  the  model  from  its  components  while  still  maintaining  a  closed-form 
representation.  In  some  cases,  expert  knowledge  can  be  used  to  assign  state- 
dependent  distributions  for  DPUT  and/or  LT. 

As  an  illustration,  suppose  the  first  row  of  10  observations  in  Table  1  can  be 
associated  with  replenishment  orders  where  a  significant  number  of  missions  were 
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canceled  due  to  weather,  creating  reduced  demand.  This  reduced  demand  is 
assumed  to  occur  on  20%  of  replenishment  orders;  thus,  demand  can  be  considered 
to  have  two  states:  regular  (with  80%  probability)  and  low  (20%  of  the  time). 

To  demonstrate  another  approach  to  finding  MTE  approximations,  the  dataset 
in  Table  1  will  be  used  in  this  example  to  first  determine  a  standard  PDF  that  best 
fits  the  empirical  data  for  each  demand  state.  In  this  case,  the  log-normal 
distribution  with  |i  =  1.03  and  =  0.31  is  selected  for  the  regular  state  and  the 
N(0.27,0.1 9)  is  chosen  for  state  2.  The  demand  in  each  state  for  a  given  lead  time 
period  is  then  a  sum  of  i.i.d.  log-normal  random  variables.  This  sum  has  no  known 
distribution,  but  approximations  for  the  PDF  of  a  sum  of  log-normal  random  variables 
exist.  Following  Cobb  et  al.  (2013),  the  Fenton-Wilkinson  approximation  (Fenton, 
1960)  is  implemented  and  MTE  distributions  are  fit  to  these  approximations  for  each 
state  and  each  possible  lead  time  value.  For  state  1  and  state  2,  these  functions  are 
denoted  by  and  respectively.  The  conditional  PDF  for  LTD  given 

L  -  Us  then  calculated  as 

fxi[L=i}M  =  0.8  ■  fx%^ij(x)  +  0.2  ■  /xfi^ijM. 

The  PDF  for  LTD  is  constructed  as  in  equation  (4).  The  new  LTD  distribution  is  bi¬ 
model,  as  shown  in  Figure  5. 


J.08 

06 

■04 

>2 

10 

Aci 

Figure  5.  Mixture  Distribution  for  LTD  With  State-dependent  Demand 

Suppose  the  state-dependent,  bi-modal  distribution  shown  in  Figure  5  is  the  correct 
PDF  for  LTD.  Using  this  distribution  as  part  of  the  total  cost  function  to  find  the 
optimal  (Q,R)  policy  results  in  a  21%  savings  when  compared  to  implementing  the 
policies  found  earlier  using  the  MTE  distribution  shown  in  Figure  4  (or  one  of  the 
other  approximations).  The  mixture  distribution  approach  still  yields  a  closed-form 
function  for  Sr  and  the  optimization  is  still  fast. 
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Coordinated  Supply  Chains 

In  this  section,  we  consider  a  two-echelon  supply  chain  as  depicted  in  Figure 
6.  A  buyer  experiencing  random  demand  places  its  orders  for  inventory  with  the 
supplier. 


Inventorx 

' 

' 

Demand 

(Random) 

Supplier 

Order 

Buyer 

Figure  6.  The  Cost  Function  for  the  Buyer  in  This  Probiem  is  as  Foiiows 

(Hadley  and  Whitin,  1963;  Johnson  and  Montgomery,  1974): 

TC,(Q.R,  F)  =  (/^,  -  K)  ■  ^  +  ^  +  /i,  ■  (0.5(2  +  R-  (7) 

Most  of  the  notation  is  the  same  as  for  the  cost  function  defined  in  equation  (5).  The 
subscript  b  has  been  added  to  the  fixed  cost  per  order,  annual  unit  holding  cost,  and 
total  cost  to  identify  this  amount  with  the  buyer.  The  subscript  s  will  similarly 
represent  the  seller.  The  quantity  V  \sa  rebate  provided  by  the  seller  to  the  buyer 
on  a  per  order  basis  as  an  incentive  for  the  buyer  to  adopt  policies  that  benefit  both 
parties  (Cobb  and  Johnson,  2014).  As  discussed  in  the  introduction,  credit  options 
and  price  discounts  have  also  been  considered  in  this  two-level  supply  chain  as 
coordination  incentives  (Chaharsooghi  and  Heydari,  2010;  Chaharsooghi  et  al., 

201 1 ;  Li  and  Liu,  2006). 

The  cost  function  for  the  supplier  in  this  problem  is: 

TC,(Q,  N,V)  =  {^  +  v)-^  +  KiN  -  1)Q.SQ.  (8) 

In  this  two-level  supply  chain  model,  the  buyer  selects  an  order  quantity  and 
reorder  point.  The  supplier  receives  orders  of  size  Qfrom  the  buyer  and  purchases 
inventory  from  its  vendors  in  a  quantity  that  is  an  integer  multiple  A/ of  the  buyer’s 
order  size. 

The  supply  chain  can  operate  in  one  of  three  modes.  First,  the  buyer  can 
select  Qc/and  F?c/ without  considering  the  effect  of  its  selection  on  the  supplier’s 
costs.  In  response,  the  supplier  selects  Nd  to  minimize  its  own  costs.  This  is 
referred  to  as  the  decentralized  mode  and  since  there  is  no  coordination,  the  rebate 
amount  is  1^=0.  Total  costs  in  the  supply  chain  are  TC^  = 

TCb{Qd,Fld,0)+TCs{Qd,Nd,0).  Second,  the  buyer  and  supplier  can  agree  on  values  for 
Qc,  Rc,  and  A/cthat  minimize  the  sum  of  the  cost  functions  in  equations  (7)  and  (8). 
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Since  the  members  cooperate  fully  and  are  centralized,  there  is  again  no 
requirement  for  the  supplier  to  provide  a  coordination  incentive  and  V=0.  Total  costs 
in  this  mode  are  denoted  by  T(y=  TCb{Qc,Rc,0)+TCs{Qc,Nc,0). 

If  the  parties  are  not  centralized  but  can  coordinate  their  policies,  the  potential 
exists  to  divide  cost  savings  of  T(T  =  TC^-TCf.  An  interval  [Vmin,Vmax]  can  be 
calculated  (Cobb  and  Johnson,  2013)  such  that  any  value  for  the  rebate  V  \n  the 
interval  reduces  the  total  costs  in  the  supply  chain  to  centralized  levels.  The 
smallest  value  of  the  rebate  the  buyer  will  accept  can  be  found  by  solving 
TCb{Qc,Rc,V)  =  TCb{Qd,Rd,0)  for  y.  This  value  is  denoted  by  t/m/n-  The  largest  value 
of  the  rebate  the  seller  will  accept  can  be  found  by  solving  TCs{Qc,Nc,  V)  = 
TCs{Qd,Nd,0)  for  V.  This  value  is  denoted  by  Vmax-  For  the  example  in  this  paper, 
we  will  assume  that  if  the  parties  agree  to  coordinate  their  policies  (and  implement 
Qc,  Rc,  and  Nc),  the  value  of  the  rebate  they  select  is  K=(  Vmin+Vmax)l2.. 

All  of  the  two-echelon  supply  chain  models  referenced  previously  assume  that 
demand  for  the  entire  lead  time  period  is  normally  distributed.  For  the  case  where 
both  Q  and  R  are  selected  to  minimize  total  costs,  Charharsooghi  and  Fleydari 
(201 0)  derive  expressions  that  state  the  optimal  value  for  Q  (in  either  the 
decentralized  or  centralized  mode)  as  a  function  of  the  optimal  value  for  R  (and  vice 
versa)  and  the  standard  normal  cumulative  density  function.  The  optimal  values  can 
be  found  by  iterating  between  these  two  expressions.  The  supplier  selects  the 
integer  value  for  A/ that  minimizes  its  costs  subject  to  the  choices  of  the  buyer. 

By  implementing  the  mixture  distribution  approach,  we  can  develop  closed- 
form  expressions  for  the  cost  functions  in  (7)  and  (8)  and  find  optimal  solutions  in  the 
same  manner  as  the  solutions  presented  earlier  in  the  paper  for  the  {Q,R)  inventory 
model.  For  illustration,  assume  Y=E{D)  ■  250  working  days  =  720,  Ks=Kb=30, 
hs=hb=4,  and  tt=5.  These  parameters  are  the  same  as  used  in  the  earlier  example 
and  the  supplier  has  the  same  cost  structure  as  the  buyer  (obviously  this  may  not 
always  be  true  in  practice). 

For  the  previous  example,  employing  the  MTE  mixture  distribution  in  Figure  4 
gives  the  same  results  in  Table  2  for  the  decentralized  case — Qd=1 10  and  Rd=27.  In 
this  mode,  the  supplier  selects  the  multiple  of  the  buyer’s  order  quantity  that 
minimizes  its  costs.  Since  rCs(1 1 0,1 ,0)=1 97  and  rCs(1 1 0,2,0)=31 6,  the  supplier 
selects  Nd=^ .  Total  supply  chain  costs  in  the  decentralized  mode  are  TC^  =  678. 

In  the  centralized  mode,  we  find  the  optimal  order  quantity  and  reorder  point 
that  minimizes  TCb{Q,R,0)+TCs{Q,N,0)  for  several  possible  values  of  N,  then  choose 
the  optimal  values  that  give  the  lowest  combined  supply  chain  cost.  Again,  using  the 
MTE  mixture  distribution  allows  the  construction  of  a  closed-form  total  cost  function, 
and  optimization  over  this  function  in  Mathematica  is  fast.  Using  the  MTE  mixture 
distribution,  we  find  that  Qc=154,  Rc=24,  and  A/c=1.  Total  supply  chain  costs  in  the 
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centralized  mode  are  TC^  =  648.  Table  3  summarizes  the  optimal  values  for  the 
decision  variables  in  each  mode  and  the  total  costs  for  each  party  and  the  supply 
chain.  The  answers  obtained  with  the  mixture  distribution  approach  are  compared 
with  those  obtained  by  using  the  solutions  shown  by  Chaharsooghi  and  Heydari 
(2010). 

Table  3.  Optimal  Solutions  and  Total  Costs  for  the  Supply  Chain  in  Three 


Modes  of  Operation 


Normal 

0 

R 

N 

V 

TCt 

TCs 

TC 

Decentralized 

108 

25 

1 

0 

483 

200 

683 

Centralized 

151 

23 

1 

0 

506 

143 

649 

Coordinated 

151 

23 

1 

8.53 

466 

183 

649 

MTE  Mixture 

0 

R 

N 

V 

TCt 

TCs 

TC 

Decentralized 

110 

27 

1 

0 

481 

197 

678 

Centralized 

154 

24 

1 

0 

507 

141 

648 

Coordinated 

154 

24 

1 

8.51 

467 

181 

648 

A  comparison  of  the  solutions  in  the  decentralized  and  centralized  models 
shows  that  the  costs  in  the  entire  supply  chain  can  be  reduced  by  TCT  =  TC^-TCf= 

30  if  the  centralized  order  quantity  and  reorder  point  are  implemented.  However, 
these  policies  increase  costs  for  the  buyer  by  507-481  =26.  By  using  the  solutions  in 
Cobb  and  Johnson  (201 3)  to  find  the  value  V  that  divides  the  cost  savings  of 
operating  in  the  centralized  mode  between  the  buyer  and  the  seller,  the  buyer  is 
adequately  compensated  for  increasing  its  order  quantity.  The  rebate  amount  for 
this  problem  is  8.51  per  order  cycle.  Both  members  experience  costs  that  are  lower 
than  in  the  decentralized  mode. 

Alternative  Approach 

This  section  introduces  an  alternative  approach  to  modeling  the  LTD 
distribution,  the  mixture  of  polynomials  (MOP)  model. 

To  illustrate  the  formation  of  the  LTD  distribution,  we  will  utilize  the  following 
example  from  McClain  and  Thomas  (1985)  that  has  also  been  used  by  Eppen  and 
Martin  (1988).  Demand  in  each  time  period  is  normally  distributed  with  mean  iiD=AQ 
and  variance  a^=30.  Lead  time  (in  periods  of  one  day)  may  take  on  the  values  7, 

12,  14,  15,  16,  and  25,  and  each  value  has  a  probability  of  1/6. 

Normal  Approximation 

Because  the  possible  values  for  LT  are  dispersed  over  the  range  from  7  to 
25,  the  distribution  for  LTD  will  be  multi-modal.  As  such,  there  is  no  one  standard 
PDF  that  is  a  good  fit.  The  typical  “textbook”  approach  to  modeling  the  LTD 
distribution  in  this  case  is  a  normal  approximation,  and  the  normal  distribution  has 
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been  used  exclusively  in  the  two-stage  supply  chain  model  under  continuous  review 
assumptions  that  will  be  presented  later  in  the  paper. 

The  normal  approximation  to  the  compound  LTD  distribution  has  a  mean  and 
variance  as  defined  in  equation  (2).  In  the  example  under  consideration,  E{L)  = 

14.83  days  and  LT  has  a  variance  of  Var{L)=29A4.  The  formulas  in  (2)  are  used  to 
determine  that  E(X)=593.33  and  Var{X)  =47 047. 2.  If  we  want  to  find  the  reorder 
point  {R)  associated  with  a  certain  service  level,  say  95%,  we  can  use  the  Excel 
function  NORM. INV(0. 95, 593.33,47047. 2^0. 5)  to  find  f?=950.  The  service  level  is 
the  probability  that  all  customer  orders  are  filled  in  given  order  cycle. 

Eppen  and  Martin  (1988)  demonstrate  that  for  this  example,  implementing 
R=950  will  actually  lead  to  very  different  service  level  than  95%.  This  is  because  the 
true  distribution  of  LTD  is  a  mixture  of  normal  distributions.  This  is  discussed  in  the 
next  section. 

Mixture  of  Normal  Distributions 

In  this  section  and  for  the  remainder  of  the  paper,  the  distribution  of  LTD  is 
denoted  by  fx-  The  distribution  of  LTD  conditional  on  a  specific  value  L  -  lior  lead 
time  is  denoted  by  fx\[L=7}-  Similarly,  the  cumulative  distribution  function  (CDF)  for 
lead  time  demand  is  denoted  by  Fx,  while  the  CDF  conditional  on  a  specific  lead 
time  L  =  Z  is  denoted  by  F;c|{l=7}- 

In  the  example  problem,  if  lead  time  is  L  =  7  days,  the  distribution  fx\{L=7}  is  a 
normal  PDF  with  mean  7-40=280  and  variance  7-30=210.  The  means  and  variances 
of  all  the  conditional  LTD  distributions  can  be  similarly  calculated.  The  marginal 
distribution  for  LTD  is  the  mixture  of  normal  distributions  calculated  as 

fx(x)  =  -ifx\[L=7](x)  +  fx\[L=12](x)  +  fx\[L=14](x)  + /a'|{L=15}  W  + 

fx\[L  =  16  }(x)  +  fxllL  =25  }W ). 

The  mixture  of  normal  distributions  for  LTD  is  shown  in  Figure  7  overlaid  on 
the  normal  approximation  with  mean  593.33  and  variance  47047.2. 
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Figure  7.  LTD  Distribution  and  Normai  Approximation 


Consider  the  reorder  point  R=950.  We  can  find  the  service  level  {SL)  associated 
with  this  reorder  point  by  evaluating  the  conditional  CDFs  Fxi{l=i]  at  950  and 
weighting  the  results  (Eppen  and  Martin,  1988).  This  is  done  as  follows 

5L(950)  =  —  (F;f|{^^7j(950)  +  F;f|{L=i2}(950)  +  F;f|{L=i4}(950)  +  F;f|{L=i5}(950) 

+  Fx|{l=i6}(950)  +  F;f|{L=25}(950)  ), 

SL(950)  =  -(1  +  1  +1  +  1  +1  +  0.034)  =  83.8%. 

6 

The  conditional  values  for  SL  given  a  certain  LT  are  calculated  using  the 
NORM.DIST  formula  in  Excel;  for  example,  the  SL  given  L=25  is 
NORM. DIST(950, 7-40, (7-30)''0. 5,1 ).  Calculation  of  the  reorder  point  associated  with 
a  desired  service  level  cannot  be  done  directly  with  the  exact  LTD  distribution,  but  a 
function  such  as  Goal  Seek  in  Excel  can  be  implemented  to  find  that  R=1014 
provides  a  95%  SL. 

Mixture  of  Polynomials  Approximation 

If  the  functional  form  of  fx  permits  closed-form  integration,  the  SL  associated 
with  a  given  reorder  point,  R,  can  be  determined  as 

R 

SL(R)  =  j  fx(x)  dx. 

0 

Since  the  functional  form  of  the  mixture  of  normal  distributions  for  the  example 
problem  cannot  be  integrated  in  this  way,  built-in  Excel  functions  for  the  normal  CDF 
were  used  to  calculate  the  service  level.  This  required  weighting  the  results  from  the 
conditional  distributions  for  each  possible  lead  time  value. 
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One  method  for  obtaining  a  closed-form  distribution  for  LTD  is  the  mixture  of 
polynomials  (MOP)  model  (Shenoy  and  West,  2011).  The  MOP  model  can  be  used 
to  approximate  PDFs  by  piecewise  polynomials  defined  on  hypercubes.  MOP 
approximations  of  standard  PDFs,  such  as  the  normal  distribution,  can  be  developed 
by  using  Lagrange  interpolating  polynomials  with  Chebyshev  points  (Shenoy,  2012). 
This  method  was  used  to  define  a  2-piece,  4th-degree  MOP  function  that 
approximates  the  standard  normal  PDF  as 

..  ^  [0.398  -  0.038Z  -  0322z^  -  0.148z^  -  0.020z^  if  -  3  <  z  <  0 
10.398  +0.038z-  0.322z2 +  0.148z^ -0.020z^  ifO  <  z  <  3. 

All  piecewise  functions  in  this  paper  are  assumed  to  equal  zero  in  undefined  regions. 
Using  this  approximation,  the  PDF  for  lead  time  demand  conditional  on  L  -  I  can  be 
determined  as 


The  MOP  function  fx  that  approximates  the  PDF  fx  for  LTD  is  determined  as 


k 


/zW  =  =  h)  ■  fx\[L=l}ix). 


i=l 


The  index  i  has  been  added  to  the  k  possible  values  for  lead  time.  This 
method  can  be  used  when  the  DPUT  distribution  is  normal,  or  at  least  in  any 
situation  where  we  are  willing  to  approximate  the  DPUT  distribution  with  a  normal 
distribution.  Notice,  this  would  be  very  different  (and  more  accurate)  than 
approximating  the  distribution  for  demand  over  the  entire  lead  time  with  a  normal 
distribution. 

For  the  example  problem,  fx  is  calculated  as 


The  MOP  approximation  to  the  LTD  distribution  is  a  relatively  compact  1 5-piece,  4'*^- 
degree  polynomial  defined  as 

fx(x) 

-20.93  +  0.333X  -  0.002x2  +  5  24  x  lO-^x^  -  5.18  x  10"^  if  236.53  <  x  <  280 

-45.50  +  0.596X  -  0.003x2  ^  g  35  ^  lO-^^x^  -  5.18  x  lO-^^x^  280.00  <  x  <  323.47 


V-261.34  +  0.996X  -  0.001x2  ^  g  Q3  ^  iQ-^x^  -  2.15  x  lO-^^x^  if  1000  <  x  <  1082.16. 

This  closed-form  function  for  the  lead  time  demand  distribution  is  easy  to 
manipulate.  It  can  be  easily  integrated  to  find  a  closed-form  function  for  the  CDF  of 
lead  time  demand  as  follows: 
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FxM 


^1059.39  -  20.93X  +  0.167x2  _  o.OOOyx^  -  1.31  X  lO-^x"*-  -  1.036  X  lO^^x^ 
2770.91  -  45.50X  +  0.298x2  _  q.OOIOx^  -  1.59  X  lO-^x"*-  -  1.036  X  lO^^x^ 
0.1667 

4621.62  -  52.22X  +  0.236x2  _  q.OOOSx^  -  5.99  X  lO-^x^  -  2.692  X  lO-^^x^ 
0.8333 


if  236.53  <  X  <  280 
280.00  <  X  <  323.47 
323.47  <  X  <  423.01 
423.01  <  X  <  480 

705.73  <  X  <  917.84 


54808  -  261.34X  +  0.498x2  _  o.0005x2  +  2.26  X  lO-^x"*-  -  4.30  X  lO^^^x^  if  1000  <  x  <  1082.16 
U  X  >  1082.16. 


Using  this  CDF  to  find  the  service  level  for  a  reorder  point  of  950  gives 

4(950)  =  5L(950)  =  83.9%. 

Evaluating  950  gives  4  at  each  possible  reorder  point  value  between  E(X)  and  the 
first  value  for  Rthat  provides  a  95%  service  level  gives  F?=1015  and  this  calculation 
requires  0.05  seconds  of  computing  time. 


In  summary,  the  LTD  distribution  can  be  modeled  using  one  normal 
distribution  as  an  approximation  over  the  entire  lead  time  period.  This  method  leads 
to  poor  results  when  calculating  the  service  level  for  a  given  lead  time  and  for  finding 
a  reorder  point  that  achieves  a  targeted  service  level.  The  actual  distribution  for  the 
example  problem  is  a  mixture  of  normal  distributions,  and  Excel  formulas  and  built-in 
functions  can  be  utilized  to  find  service  levels  and  reorder  points,  albeit  indirectly. 
The  MOP  model  offers  an  alternative  to  constructing  a  closed-form  LTD  distribution 
that  can  be  directly  integrated  and  evaluated  to  find  a  CDF  for  lead  time  demand, 
service  levels,  and  reorder  points.  As  discussed  in  the  remainder  of  the  paper,  this 
distribution  can  be  utilized  to  find  optimal  inventory  policies  in  a  two-level  supply 
chain  under  uncertain  demand  and  continuous  review  assumptions. 


To  implement  the  MOP  mixture  distribution  approach  to  find  an  optimal  order 


quantity/reorder  point  combination,  we  first  develop  a  closed-form  expression  for  the 


expected  shortage  per  cycle  in  (6)  using  the  previously  defined  PDF  4.  This 

function  is  an  8-piece,  G'^^-degree  polynomial  defined  as 

Sr(R) 


f-3.16  X  10“^  +  32362AR  -  138.1R^  +  0.314i?^ 
-0.0004i?4  +  2.75  X  -  7.81  x 

-4.11  X  10“^  +  40257.0i?  -  164AR^  +  0.358i?^ 
-0.0004«^  +  2.87  X  IQ-^R^  -  7.81  x 


if 593.33  <R  <600 
if  600  <R<  621.48 


-9.57  X  10-^  +  54807.6i?  -  130.7i?2  +  0.166i?^ 

V  -O.OOOli?^  +  4.52  X  10-^i?^  -  7.16  x  10-^^  R^  if  1000  <R<  1082.16. 


Decentralized  Solution 

This  function  for  5^  shown  above  can  be  substituted  into  equation  (5)  to 
create  a  piecewise  cost  function  for  the  buyer.  In  this  example,  we  will  assume 
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Kh  -  SO,hi,  =  5,  and  tt  =  6.  Expected  annual  demand  is  based  on  150  working  days 
and  equals  Y  =  150  •  \iD  =  150-40  =  6000.  This  cost  function  is  displayed  as  a 
function  of  Q  for  three  values  of  R  in  Figure  8.  By  inspection,  we  can  see  that  the 
optimal  order  quantity  is  lower  for  smaller  values  of  R.  In  other  words,  we  can 
better  control  costs  by  simultaneously  selecting  the  order  quantity  and  reorder  point. 


Figure  8.  Buyers’  Cost  as  a  Function  of  Order  Quantity  for  Three  Vaiues  of 

the  Reorder  Point 

Optimization  over  the  cost  function  developed  using  the  MOP  distribution  for  LTD  is 
fast.  Notice  that  the  function  for  expected  shortage  per  cycle  is  a  MOP.  When  this 
expression  is  inserted  in  the  cost  function  in  (5),  the  result  is  a  function  with 
polynomial  terms  and  some  terms  with  Q  in  the  denominator.  The  example  here 
was  solved  using  Mathematica  9.0  by  using  the  ArgMin  function.  The  resulting 
solutions  are  0^=364  and  Rd=1014  with  TCb{Qd, Rd, 0)=3924.  The  supplier’s  best 
response  is  to  set  Nd=^  and  incur  costs  of  TCs{Qd,Nd,0)=2472,  and  total  costs  in  the 
supply  chain  are  rc°'=6396.  The  computing  time  expended  is  less  than  one  second. 

An  iterative  approach  (Hadley  and  Whitin,  1961)  in  combination  with 
numerical  integration  was  implemented  to  find  the  solutions  using  the  normal 
approximation  to  the  LTD  distribution  using  the  partial  solution  provided  by 
Chaharsooghi  and  Heydari  (2010).  The  solutions  are  =  447  and  =  925.  If 
these  solutions  are  inserted  in  the  “actual”  cost  function  (the  one  developed  with  the 
MOP  distribution  for  LTD),  the  result  is  =  4454.  Using  the  MOP 

mixture  distribution  yields  an  improvement  in  costs  of  4454  -  3924  =  530  or  12%. 

Centralized  Solution 

The  closed-form  function  Spfor  expected  shortage  per  cycle  developed  using 
the  MOP  distribution  for  LTD  can  also  be  used  to  derive  a  cost  function  for  the  entire 
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supply  chain  in  the  centralized  case.  In  this  example  we  assume,  Ks=^  50  and 
hs=12.5.  This  function  TC^[Q,R,N)  is  used  to  find  the  optimal  combination  (Qc,  for 
several  possible  values  of  the  supplier’s  decision  variable  N.  The  value  of  N 
producing  the  lowest  total  cost  once  the  corresponding  optimal  values  for  order 
quantity  and  reorder  point  are  selected  is  deemed  the  best  supplier  policy.  Typically, 
solving  for  the  optimal  (Qc,  R^  with  A/=1  then  checking  to  see  if  N=2  or  A/=3 
produces  a  better  solution  is  adequate. 

The  best  order  quantity  in  the  centralized  model  for  a  given  reorder  point  is 
higher  than  the  optimal  order  quantity  in  the  decentralized  case.  This  is  illustrated  in 
Figure  9,  where  the  total  costs  are  graphed  as  a  function  of  Q  for  the  decentralized 
and  centralized  cases  assuming  a  reorder  point  of  R=^000.  Visually,  the  centralized 
cost  function  appears  to  reach  a  minimum  at  a  larger  value  of  Q. 


Figure  9.  Decentralized  and  Centralized  Costs  as  a  Function  of  Order 

Quantity  for  a  Reorder  Point  of  1000. 

The  closed-form  centralized  cost  function  can  again  be  easily  utilized  to  find  the 
optimal  policy  of  {Qc,Rc,Nc)  =  (718,993,1).  The  costs  for  the  parties  at  the  optimal 
solutions  are  as  follows:  TCb{Qc,Nc,0)=4333-,  rCs(Qc,A/c,  V)=1254;  TCf =5587.  The 
solution  again  takes  around  one  second  of  computing  time  to  obtain. 

The  buyer  incurs  higher  costs  by  4333-3924=409  in  the  decentralized  mode 
as  compared  to  the  centralized  mode,  where  the  supplier’s  costs  are  reduced  by 
2472-1 254=1 21 8.  Total  costs  in  the  supply  chain  are  lower  than  in  the  decentralized 
mode  by  6396-5587=809. 

The  corresponding  centralized  solutions  found  using  the  normal 
approximation  are  -  802  and  Rfj  -  857.  If  these  solutions  are  inserted  in  the 
“actual”  cost  function  for  the  supply  chain  (the  one  developed  with  the  MOP 
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distribution  for  LTD),  the  result  is  TC^  =  5891 .  Using  the  MOP  mixture  distribution 
yields  an  improvement  in  costs  of  5891  -  5587  =  304  or  5%  in  the  centralized  mode. 

Coordinated  Solution 

While  the  buyer  would  prefer  that  the  supply  chain  operate  in  decentralized 
mode  and  the  supplier  wants  a  centralized  solution,  both  parties  can  potentially 
compromise  and  coordinate  to  divide  the  centralized  costs  savings.  The  closed-form 
cost  functions  developed  using  the  MOP  method  again  provide  an  approach  to 
determine  a  supply  chain  coordination  mechanism  to  make  this  work. 

The  buyer  will  accept  a  per  order  rebate  as  low  as  Vmin,  which  can  be  found 
by  solving  70^(718,993,  V)  =  3924,  or  4333-8. 356v=3924.  The  solution  is  Vmin=^9. 
The  supplier  will  accept  a  per  order  rebate  as  high  as  Vmax,  which  can  be  found  by 
solving  7Cs(71 8,993,  =  2472,  or  1 254-1-8. 356v=2472.  The  solution  is  V'max=146. 

In  this  example,  at  the  centralized  optimal  order  quantity,  there  are 
V7Qc=6000/718=8.356  order  cycles  per  year,  so  the  minimum  incentive  entails 
rebates  of  8.356-49=409  and  the  maximum  incentive  entails  rebates  of 
8.356-146=1218.  One  solution  is  to  implement  V={Vmin+Vmax)/2=97.5  and  require 
the  supplier  to  provide  815  in  rebates  to  the  buyer.  This  brings  the  buyer’s  total 
costs  to  3518,  the  supplier’s  total  costs  to  2069,  and  supply  chain  costs  to  5587, 
which  is  the  centralized  level. 

Empirical  MOP  Distributions 

In  previous  sections  of  the  report,  we  have  seen  the  MTE  approach 
implemented  with  empirical  data,  and  the  MOP  approach  implemented  when  the 
underlying  LT  and  DPUT  distributions  were  discrete  and  represented  by  a  standard 
continuous  PDF,  respectively.  This  section  will  illustrate  an  approach  to  estimating 
an  MOP  function  to  approximate  the  LTD  distribution  when  empirical  data  is 
available. 

We  suppose  a  modest  amount  of  historical  data  is  available  for  daily  demand 
and  lead  times.  In  this  example,  we  use  a  dataset  Doi  N=  500  observations  for 
daily  demand  with  sample  mean  39.66  and  sample  variance  30.64.  These  values 
are  a  random  sample  from  the  A/(40,30)  distribution.  Fifty  observations  {Nl=50)  of 
historical  lead  time  values  are  available  in  dataset  DLwith  sample  mean  14.94  and 
sample  variance  37.04.  There  are  values  in  the  dataset  for  each  possible  lead  time 
value.  Empirical  histograms  are  displayed  for  this  sample  data  in  Figure  10. 
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Figure  10.  Histograms  for  the  Daiiy  Demand  and  Lead  Time  Data. 

The  datasets  depicted  in  Figure  10  will  be  used  to  develop  LTD  distributions 
given  each  possible  empirical  lead  time  value.  We  assume  that  the  observations  of 
daily  demand  are  i.i.d.  Since  this  is  the  case,  we  use  the  dataset  oi  N=  500  values 
to  create  six  smaller  datasets  for  daily  demand  given  the  possible  values  for  L  For 
example,  the  first  seven  values  for  daily  demand  are  40,  37,  47,  47,  45,  34,  and  31 . 
These  are  summed  to  281  to  determine  the  first  sample  value  in  the  dataset  for  LTD 
demand  given  a  lead  time  of  7  days.  The  next  seven  consecutive  values  in  the 
dataset  sum  to  263,  so  this  is  the  second  value  in  the  L=7  dataset,  and  so  on.  This 
smaller  dataset  has  71  observations. 


B-Spline  Estimation  of  MOPs 

Lopez-Cruz  et  al.  (2014)  suggest  using  a  linear  combination  of  B-spline 
functions  to  construct  MOP  approximations  from  datasets  where  the  parametric  form 
of  the  underlying  probability  distribution  is  unknown.  B-spline  functions  are 
piecewise  polynomial  functions  defined  by  the  number  of  control  points,  n-i-1 ,  and  the 
degree  of  the  polynomial,  d.  The  control  points  define  a  knot  vector  t  =  {to,U,t2,  ■■■, 
fn}. 


B-spline  functions  (Zong  and  Lam,  1998)  have  two  definitions,  one  when  d=^ 
and  another  when  d>^ .  When  d=^ ,  the  functions  are  defined  as 


tj<X  <  tj+^ 

otherwise  . 


For  d>^ ,  the  functions  are  calculated  as 


X  —  ti 


d'j+k-l  / 


^j+k  ^ 
S-j+k  ~  ^j+1/ 


Bj  +  i,k-l(.^y 


The  control  points  are  indexed  by  y=0,...,nand  the  degree  of  the  functions  are 
indexed  by  k=^,...,d.  For  this  example,  we  assume  k  is  the  smallest  value  in  the 
dataset,  tn  is  the  largest  value  in  the  dataset,  and  that  the  intervals  between  all  of  the 
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knots  are  the  same  distance.  The  resulting  functions  Bj^k  are  referred  to  as  uniform 
B-splines.  The  B-spline  functions  are  used  to  form  a  n-piece  MOP  density  function 

m 

fx\L=i(.x)  =  ^  (x)  (9) 

i=l 


by  selecting  mixing  coefficients  a,,  where  m=n+d-^.  Thus,  the  PDF  for  LTD 

given  a  lead  time  /.=/will  be  a  mixture  of  the  m  B-splines  of  order  d. 


Suppose  a  dataset  D  =  {x^,  ...,x„^ }  of  observations  of  LTD,  X,  given  a  specific 


lead  time  L  -  r\s  available.  Zong  (2006)  suggests  using  the  following  iterative 
formula  for  determining  the  maximum  likelihood  estimators,  {a^, ... ,  for  the 
mixing  coefficients  in  (9): 


d 

Ni  ■  (ti  -  tj_d) 


I 

X£D 


fx\L=i  [x]  . 


(10) 


Beginning  with  equivalent  values  for  each  a,,  the  expression  in  (10)  is  used 
iteratively  for /=1,...,m  until  f  where  is  the  log-likelihood 

of  D  given  fx\L=l  at  iteration  q  in  the  optimization  process.  Using  £=10  ® 

appears  to  be  adequate  for  most  applications  (Lopez-Cruz  et  al.,  2014). 

The  goal  is  to  develop  PDFs  fx\L=i  for  LTD  given  each  possible  value  for  L 
that  are  reasonably  accurate;  however,  we  would  like  the  number  of  pieces  and  the 
degree  of  the  polynomial  functions  comprising  the  MOP  densities  to  be  as  small  as 
possible  to  avoid  overfitting  and  speed  up  computation  of  optimal  inventory  policies. 
Thus,  we  will  consider  several  possible  values  for  c/and  nfor  each  PDF  and  select 
the  approximation  that  maximizes  the  Bayesian  information  criterion  (BIG)  calculated 
as 

BIC{fx\L^l(ix),D)  =  L(£)|/;,|,.j(x))  (11) 


The  second  term  in  the  BIG  expression  is  a  penalty  for  adding  parameters  to 
the  model.  The  approach  we  will  take  is  to  find  the  values  of  c/and  n  that  maximize 
the  BIG  score  for  the  PDF  of  LTD  given  the  most  likely  value  L=l,  then  use  those 
parameters  to  estimate  each  of  the  conditional  PDFs  for  LTD  given  each  possible 
lead  time  value.  In  practice,  once  we  settle  on  good  values  for  c/and  n,  this  step 
could  be  avoided.  Alternatively,  we  have  found  that  d=3  and  n=3  seem  to  be 
adequate  for  many  problems. 

Notice  from  Figure  10  that  L=7  is  the  lead  time  value  that  occurs  most 
frequently  in  the  empirical  data.  Thus,  we  will  begin  by  constructing  the  PDF  fx\L=7 
for  LTD  given  L=7  with  an  MOP  function  constructed  from  B-splines  with  the  values 
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of  c/and  n  that  maximize  the  BIG  score.  We  tested  all  possible  combinations  of 
parameter  values  where  c/and  n  are  in  the  set  {2, 3, 4, 5, 6}. 

The  best  value  for  the  BIG  score  is  -290.4  and  is  achieved  for  n=2  and  c/=3. 
The  PDF  is  shown  graphically  in  the  left  panel  of  Figure  1 1  overlaid  on  the 
A/(7-40,7-30)  distribution.  Recall  that  the  MOP  is  not  fit  to  the  normal  PDF,  but 
rather  a  small  sample  of  data  generated  from  the  normal  PDF.  The  four  B-splines 
used  to  construct  the  MOP  function  are  shown  in  the  right  panel  of  Figure  1 1 .  The 
mixing  coefficients  determined  via  37  iterations  of  equation  (10)  are  ai=0.04, 
02=0.07,  03=0.83,  and  04=0. 06. 


Figure  11.  PDF  for  LTD  Given  a  Lead  Time  of  7  days  (ieft)  and  B-spiines 
Used  to  Construct  the  MOP  Function  (right) 


The  MOP  fitting  process  is  repeated  for  /=1 2,1 4,1 5,1 6,25  to  find  each 
conditional  PDF  fx\L=i-  The  mixture  distribution  for  LTD  is  then  calculated  (Gobb, 
2013)  as 

6 

fx(x)  =  Y^P{L  = 

'  f  X\L=l^'^^ 

The  superscript  {k)  has  been  added  as  an  index  for  the  number  of  possible  LT 
values.  Here,  k=^  corresponds  with  L  =  =  7,k  =  2  corresponds  with  L  =  - 

12,  and  so  on.  The  resulting  distribution  fx  is  displayed  in  Figure  12  overlaid  on  the 
actual  distribution  for  LTD.  Again,  fx  was  created  from  the  sample  data  without 
knowledge  of  the  underlying  LTD  distribution.  The  function  fx  is  relatively  compact- 
-it  contains  1 5  pieces  and  is  a  2nd  degree  polynomial.  Since  the  class  of  MOP 
functions  is  closed  under  addition  and  multiplication,  the  resulting  mixture  distribution 
is  also  an  MOP  function. 
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Figure  12.  MOP  Mixture  Distribution  for  LTD  Overiaid  on  the  Actuai  LTD 

Distribution. 

The  MOP  distribution  for  LTD  is  utilized  to  find  a  closed-form  expression  for 
expected  shortage  per  cycle  as 

^max 

=  J  (x-R)-/xM  dx. 

R 

The  value  Xmax  is  the  largest  value  for  which  is  defined  as  non-zero.  In  the 
previous  example,  5^  is  a  9-piece,  5th  degree  polynomial.  Again,  both  {x-R)  and  fx 
are  MOP  functions,  so  the  resulting  expression  for  5^  is  an  MOP  function  because 
this  class  of  functions  is  closed  under  multiplication  and  integration.  Since  has  a 
closed-form,  the  cost  function  TCb  can  be  derived  and  used  to  find  optimal  inventory 
policies.  The  function  TCb  is  a  9-piece  function  with  polynomial  terms  in  Q  and  R 
and  some  terms  that  have  a  polynomial  numerator  and  a  Qterm  in  the  denominator. 

An  outline  of  the  process  for  using  empirical  data  on  DPUT  and  LT  to 
construct  the  LTD  distribution  fx  and  buyer's  cost  function  is  as  follows: 

1 .  Collect  datasets  of  observations  for  DPUT  and  LT.  Each  of  the 
possible  observations  for  LT  form  the  set  Ql  of  values  ..., 

2.  Create  datasets  for  LTD  given  each  possible  LT  value  by  summing 
consecutive  values  as  many  times  as  possible  from  the  DPUT  dataset. 

3.  For  the  most  likely  LT  value  in  the  discrete  empirical  distribution,  use 
the  corresponding  LTD  dataset  to  calculate  the  MOP  distribution 
fx\L=iW  that  maximizes  the  BIC  score  in  (1 1 )  by  testing  different 
values  of  c/and  nand  using  equation  (10)  iteratively  until  convergence. 

4.  Using  the  values  of  c/and  nthat  maximize  BIC  in  the  previous  step  for 
L=  create  MOP  distributions  fx  for  each  additional  value  in  Q/.. 
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(12) 


5.  Calculate  the  LTD  distribution  as 

K 

k  =  l 

6.  Calculate  the  expected  shortage  per  cycle  function  and  insert  this  in 
the  buyer's  cost  function  to  create  TCb. 

Optimal  Policies  and  Results 

This  example  again  assumes  Kb=50,  Ks=^50,  lib=5,  hs=^2.5,  and  tt=6.  We 
assume  250  working  days  per  year  so  Y  =  250-39.66  «  10000.  The  results  obtained 
for  optimal  inventory  policies  using  the  MOP  approximation  developed  from 
empirical  data  using  the  B-spline  estimation  method  are  again  compared  to  the 
solution  developed  by  Chaharsooghi  and  Heydari  (2011)  (the  CH  solution)  that 
assumes  a  normal  distribution  for  LTD.  The  total  costs  are  obtained  by  simulating 
the  expected  total  costs  on  10,000  simulation  trials  using  the  actual  normal 
distributions  for  daily  demand  and  the  discrete  distribution  for  lead  time.  Clearly,  the 
cost  savings  are  dependent  on  the  parameters. 


Table  4.  Optimal  Solutions  and  Total  Costs  for  the  Supply  Chain 


Decentralized 

0 

R 

N 

TCb 

%  Dec. 

TC 

CH  (Normal) 

558 

999 

1 

4531 

- 

683 

MOP  (Mixture  Dist.) 

455 

1019 

1 

4455 

2.0% 

649 

Coordinated 

Q 

R 

N 

TCf 

%  Dec 

TC 

CH  (Normal) 

1012 

926 

1 

6872 

- 

678 

MOP  (Mixture  Dist.) 

909 

1004 

1 

6628 

3.7% 

648 

Conclusions 


This  paper  serves  as  an  introduction  to  using  a  mixture  distribution  approach 
to  modeling  the  probability  density  function  for  lead  time  demand  in  problems  where 
a  continuous  review  inventory  system  is  implemented.  First,  construction  of  the  lead 
time  distribution  was  illustrated.  This  distribution  was  then  utilized  to  determine 
optimal  order  policies  in  cases  where  a  buyer  makes  its  decisions  alone,  and  then 
when  members  of  a  two-level  supply  chain  coordinate  their  actions. 

Several  approaches  to  modeling  the  LTD  distribution  were  illustrated.  We 
first  considered  using  the  normal  CDF  for  each  possible  lead  time  value  and 
weighting  the  results  with  discrete  probabilities.  This  method  allows  calculation  of 
service  level  probabilities,  but  does  not  provide  a  closed-form  approximation  to  the 
LTD  distribution  that  can  be  used  for  determining  optimal  inventory  policies.  Next,  a 
mixture  of  truncated  exponentials  approximation  was  utilized  to  model  the  LTD 
distribution.  The  distributions  for  each  possible  lead  time  value  were  estimated,  then 
the  results  were  weighted  to  create  a  mixture  distribution  for  LTD.  Finally,  mixture  of 
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polynomial  approximations  were  considered.  The  use  of  this  model  was 
demonstrated  as  an  approximation  to  both  standard  probability  density  functions  and 
as  a  fit  to  actual  historical  data.  Use  of  either  approximation  technique  allows 
determination  of  optimal  inventory  policies  and  provides  significant  cost  savings  as 
compared  to  a  solution  where  demand  over  the  entire  lead  time  period  is  assumed 
to  be  normal. 
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